Fragile-to-strong transition and polyamorphism in the energy landscape of liquid silica 
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We study the static and dynamic properties of liquid silica over a wide range of temperature T and 
density p using computer simulations. The results reveal a change in the potential energy landscape 
as T decreases that underlies a transition from a fragile liquid at high T to a strong liquid at low T. 
We also show that a specific heat anomaly is associated with this fragile-to-strong transition, and 
suggest that this anomaly is related to the polyamorphic behaviour of amorphous solid silica. 



Liquid silica is the archetypal glass former, and com- 
pounds based on silica are ubiquitous as natural and 
man-made amorphous materials. Liquid silica is also the 
extreme case of a "strong liquid," one having a variation 
of viscosity r\ with temperature T that closely follows the 
Arrhenius law log r/ <~ 1 /T as the liquid is cooled toward 
its glass transition temperature T g |l} |2j. In contrast, 
most liquids are to some degree "fragile," exhibiting sig- 
nificantly faster-than- Arrhenius increases of r\ as T — > T g . 
Recent studies focusing on the properties of the potential 
energy hypersurface (PES) — or "energy landscape" — of 
the liquid have demonstrated the controlling influence of 
the PES on transport properties as T -»■ T g ||, g, |, 
However, these studies have only addressed fragile liq- 
uids, and the origin of strong liquid behaviour in terms 
of the energy landscape has not yet been resolved. It is 
this question that we address here. 

In a molecular dynamics computer simulation of an 
equilibrium liquid, the diffusion coefficient D is readily 
evaluated from the particle trajectories. Like 77, D is a 
characteristic transport property whose deviations from 
the Arrhenius law serve to classify a liquid as strong or 
fragile. The theory of Adam and Gibbs (AG) [|| states 
that D is related to the configurational entropy, S c via, 

D = D exp(-A/TS D ), (1) 

where the parameters Do and A are commonly assumed 
to be T-independent. S c quantifies the number of distinct 
configurational states explored by the liquid in equilib- 
rium. These states have been proposed to correspond to 
the "basins" of the PES sampled by the liquid [§, @. A 
basin is the set of points in phase space representing con- 
figurations having the same local minimum. The local 
minimum configuration is termed an inherent structure 
(IS), and is identified in simulation by a steepest descent 
minimization of the potential energy. 

Following the IS thermodynamic formalism of Still- 
inger and Weber Q, the internal energy of the liquid 
can be expressed as E = ejs + Eharm + E an h, where 
e/s is the average IS energy and the last two terms are 
the average contributions to E due to thermal excita- 
tions about the IS. Eharm is the average harmonic con- 
tribution determined from a quadratic approximation to 
E around each IS minimum, while E an h is the remain- 



ing, necessarily anharmonic contribution. The harmonic 
and anharmonic potentials characterize the shape of the 
basin. 

If the shape of the basins does not depend on e/s (a 
condition satisfied at constant density p in the present 
study; see Methods), then S c can be calculated along 
an isochore by integrating the T dependence of e/s at 
constant volume V [Q: 

where S® is the value of S c at a reference T = Tq (see 
Methods) . Eq. || highlights that the T dependence of S c 
arises solely from changes in e/s |,|. Evaluation of S c 
when the basin shape depends on e/s is also feasible Iq, 




For a strong liquid that satisfies Eq. [I], the T depen- 
dence of S c must approach a constant to recover Arrhe- 
nius behaviour. From Eq. ||, if S c is constant, then so 
must be the variation of e/s with T. This behaviour 
would be qualitatively different from that found in sim- 
ulations of fragile liquids. For example, recent stud- 
ies of a binary Lennard- Jones liquid have shown that 
e/s oc —1/T, a dependence that results from a Gaus- 
sian distribution of IS energies §, §. When the IS en- 
ergy distribution is Gaussian, fragility has been shown 
to depend on the total number of IS states, the width 
of the Gaussian, and the dependence of the basin shape 
on e/s H However, it is not known if the PES of a 
strong liquid is similarly characterized by a Gaussian dis- 
tribution of IS energies but with parameters that differ 
from the fragile liquid case, or if the PES is qualitatively 
different from that of a fragile liquid. 

In addition, some analyses of experimental data for 
silica suggest that the liquid may be fragile at very high 
T Recent simulations @ of the BKS 0] model 

of silica between T = 2750 and 6000 K have also shown 
that at the onset of slow dynamics, as reflected for exam- 
ple by the presence of two-step relaxation in structural 
autocorrelation functions, BKS silica is a fragile liquid. 
At about 3300 K BKS silica transforms to a strong liq- 
uid and the T-dependence of all characteristic relaxation 
times becomes Arrhenius [ H . The relationship of the 



2 



-1.85 



-1.86 



-1.87 



-1.88 



Op=3.01 g/cm 
• p=2.36 g/cm ; 




0.1 0.2 0.3 0.4 
1000/T 



1000 2000 3000 4000 5000 6000 7000 8000 
T(K) 

FIG. 1: eis as a function of T along two isochores. At 
T — we show the energy Eq of the crystalline system at 
p — 2.36 g/cm 3 (filled square) and p = 3.01 g/cm 3 (open 
square). Eq is found by calculating the V dependence of the 
potential energy at T = of three crystal polymorphs of silica 
(stishovite, coesite and quartz), and then using the common 
tangent construction to determine the potential energy of the 
heterophase of coexisting crystals that would be the ground 
state at the required bulk value of p. Inset: eis for the same 
isochores as in the main panel, plotted as a function of 1/T. 
Only the high density data show a clear 1/T behaviour at low 
T. 



PES to such a "fragile-to-strong" transition is not yet 
known. 

To address these questions, we conduct extensive equi- 
librium simulations of BKS silica over a large range of 
T and p to examine the fragile-to-strong transition, and 
to identify the energy landscape behaviour that under- 
lies it. The results clarify our understanding not only of 
the origins of silica's status as a strong liquid, but also of 
the dynamical behaviour of a wider class of liquids that 
are to some degree silica-like, most notably water and 
silicon |L8j, but also other systems such as BeF2 Jl9| . 
Indeed, the concept of a "fragile-to-strong" transition 
was first proposed for the case of deeply supercooled wa- 
ter @. 

Fig. [l] shows the T dependence of ejg along two iso- 
chores. The shape of the higher p isochore is simi- 
lar to that found for fragile liquids. However, at the 
lower p — which is close to the experimental p at ambient 
pressure — eis exhibits an inflection, passing from con- 
cave downward at high T to concave upward at low T. 
Fig. [l] also shows the potential energy Eq at T = of 
the corresponding equilibrium crystalline system for the 
same p. Since the energy of the lowest IS cannot be less 
than Eq, Eq sets a lower bound on eis- The value of 
Eq relative to the measured eis curves confirms that an 
inflection occurs. 

The inset of Fig. [I] shows that the relation eis c« —1/T, 
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FIG. 2: S and its component contributions as a function of 
T at p = 3.01 g/cm 3 (left panel) and p = 2.36 g/cm 3 (centre 
panel) . S c as a function of T along two isochores (right panel) . 
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FIG. 3: log D versus 1/TS C along two isochores. Inset: Arrhc- 
nius plot of 1/D along the ambient pressure isobar, found by 
interpolating our isochoric data; and along the p — 2.36 g/cm 3 
isochore. Note the non-Arrhenius to Arrhenius crossover ob- 
served both at constant p and at constant P. All D values 
reported here are for Si atoms. 



the hallmark of a Gaussian distribution of IS energies, is 
not obeyed along our lower p isochore. The breakdown 
of this relation does not arise from changes in the shapes 
of the basins sampled at different T. Indeed, the distri- 
bution of curvature of the PES around the IS (the vibra- 
tional density of states) is found to be independent of the 
basin depth for all the isochores we study. Hence, within 
the harmonic approximation, all basins of the same p are 
characterized by the same vibrational entropy. 

Like eis, S c (Eq. ||) also exhibits an inflection along our 
lower p isochore (Fig. 0). In the same range of T, we find 
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(Fig. ||) that along each isochore D satisfies the AG rela- 
tion within numerical error. That is, despite the change 
in the nature of the T dependence of S c , the transport 
properties of the liquid adjust (from non-Arrhenius to Ar- 
rhenius) so as to maintain the AG relation, a finding that 
demonstrates of the robustness of Eq. |] ^2| . The low 
p data thus reveal the signature in the energy landscape 
of a fragile-to-strong transition. The rapid variation of 
landscape properties at high T corresponds to a fragile 
regime. As T decreases, the inflections of eis and S c 
mark the onset of a regime in which the rate of change of 
these quantities decreases, consistent with the system's 
approach to the strong liquid limit. 

We note that for the present model, S c is much smaller 
than the vibrational entropy (Sharm + S an h) and of the 
same order as S an h- Moreover, we find that for BKS 
silica the crystalline vibrational entropy differs signifi- 
cantly from the liquid vibrational entropy, as also found 
by Richet Q in his extensive analysis of silicate melts. 
This confirms that for silica, the identification of the ex- 
cess entropy (liquid entropy minus crystal entropy) with 
S c may not be adequate jg^, This highlights the 

value of finding S c via the "all-liquid route" j25) imple- 
mented here. 

The influence of the energy landscape is sufficiently 
prominent to appear in the total thermodynamic prop- 
erties |26|]. The isochoric specific heat Cy can be writ- 
ten as Cy = C v s + 3i? + C v nh , where each term is the 
derivative with respect to T of the corresponding contri- 
bution to E, and R is the gas constant. The inflection in 
the T dependence of eis corresponds to a maximum in 
C{? (Fig. H) that is the origin of a Cy anomaly, in the 
form of a peak, in the interval of T corresponding to the 
fragile-to-strong transition. An analogous Cy anomaly 
has recently been predicted to occur in the silica-like liq- 
uid BcF2 fl9|] , and in theoretical models designed to re- 
produce a fragile-to-strong transition [ p7fl . High T ex- 
periments that test for this anomaly, though challenging, 
can thus directly seek the thermodynamic signature of 
the fragile-to-strong transition in these systems. 

The peak we find in Cy occurs at T near the tem- 
perature of maximum density of BKS silica, and in the 
vicinity of a maximum of the isothermal compressibil- 
ity Kt predicted for this model pq |. Thermodynamic 
anomalies, such as peaks in Cy and Kt, have been as- 
sociated with the physics of liquid-liquid transitions and 
with polyamorphism in glasses, i.e. the abrupt pressure- 
induced transition of a low-density glass to a high-density 
glass [29 . Polyamorphism is observed experimentally 
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during compression of silica glass, and the Cy maximum 
found here may be the thermal anomaly corresponding 
to polyamorphism's mechanical anomaly. Along different 
isochores we find that the T at which the Cy maximum 
occurs decreases with increasing p, as would be expected 
for an anomaly related to polyamorphism in silica. 

These observations together suggest that the present 
fragile-to-strong transition is the dynamical transition 
corresponding to the thermodynamic crossover in the liq- 
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FIG. 4: Upper panels: Isochores of the liquid potential energy 
U. Lower panels: Cv — (3/2)R (solid lines) and C v s (dashed 
lines) as a function of T. In the upper panels are shown 
lines of slope §-R whose values at T — are the potential 
energies of the corresponding crystalline systems, calculated 
as described in Fig. |l| these lines are lower bounds for U of the 
liquid. In the left-most lower panel are shown (filled squares) 
estimates of the configurational part of Cv recently evaluated 
by Scheidler, et al. pi| , calculated from the time dependence 
of the temperature fluctuations. 



uid that presages polyamorphism p7| . More generally, 
our results provide a basis for considering all strong liq- 
uids as candidates for polyamorphism: a strong liquid 
arises via a fragile-to-strong transition, associated with 
which may be thermodynamic anomalies that are the 
liquid-state precursors to polyamorphism. 

Appendix: Methods 

Computer simulations: Results are based on molecular 
dynamics simulations of BKS silica. All data are ob- 
tained from systems of N = 1332 atoms (888 O and 444 
Si atoms), except for the p = 2.36 g/cm 3 isochore, where 
we employ 999 atoms to facilitate equilibration at low 
T. At p = 2.36 g/cm 3 8 independent runs for each state 
point are performed. All data reported here are for liq- 
uid states in equilibrium. Equilibration is confirmed by 
ensuring that runs are longer than the slowest relaxation 
time in the system as evaluated from the collective (co- 
herent) density-density correlation function. The lowest 
T runs exceed 80 ns. Simulations are carried out in the 
constant- (A, V, E) ensemble, and long-range electrostatic 
interactions are accounted for by Ewald summation. We 
evaluate eis by conducting conjugate gradient minimiza- 
tions of up to 1000 equilibrium liquid configurations and 
averaging the results. 
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Evaluation of the total entropy: For a given p, S at 
T = T = 4000 K is calculated using thermodynamic 
integration: We first find the free energy difference be- 
tween an ideal gas and a binary mixture Lennard-Jones 
(LJ) system at the chosen p by integrating the excess 
pressure along an isotherm from the high V limit where 
the system behaves as an ideal gas. We then carry out 
a set of simulations at constant V and T that continu- 
ously convert the LJ system to BKS silica, by employing 
a hybrid potential <j> = \4>bks + (1 — ^)4>l.j An 
appropriate thermodynamic integration from A = to 1 
yields the free energy of BKS silica, from which S at To 
and the chosen p is calculated. S at other T is found by 
further thermodynamic integration at constant V . 

Evaluation of the vibrational entropy: We evaluate Sharm 
from the spectrum of eigenfrequencies v (i.e. the vi- 
brational density of states) calculated from the IS's at 
each T, using S harm = fc^tJl ~ log(ft^/fcT)], where 
k and h are Boltzmann's and Planck's constants, re- 
spectively. To obtain S an h we use E an h. We evaluate 
Eanh = E-E harm -eis and then fit E anh with a polyno- 
mial constrained to be zero and have zero slope at T = 0. 
When the shape of the basins does not depend on ejg, 
Sanh may be calculated by thermodynamic integration 



using the fitted form of E an h from T = to the desired 
T. In terms of the above quantities, the integration con- 
stant in Eq. |is thus 5° = S(T ) - S anh (T ) - S harm {T ). 

Isochoric invariance of basin shape: Our assumption that 
the shape of the basins is independent of ejs along an 
isochore is based on two observations: (i) We find that 
the vibrational density of states (the v spectrum) is in- 
dependent of e/s along an isochore. (ii) The anharmonic 
energy of IS configurations heated from T = to a chosen 
T follows the Eanh = E — Eharm — e/s curve calculated 
from equilibrium simulations. This is only possible if the 
anharmonic character of the basins is also independent 
of e/s- 
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